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ABSTRACT 

The secular approximation for the evolution of hierarchical triple configurations has proven to be 
very useful in many astrophysical contexts, from planetary to triple-star systems. In this approxima- 
tion the orbits may change shape and orientation but the scmimajor axes are constant. For example, 
for highly inclined systems, the Kozai-Lidov mechanism can produce large-amplitude oscillations of 
the eccentricities. Here we re-derive the secular evolution equations including both quadrupole and oc- 
tupole orders using Hamiltonian perturbation theory. Our new derivation corrects an error in previous 
treatments of the problem. Already to quadrupole order, our results disagree with the previous "stan- 
dard" treatment; they agree only in the test-particle limit where one of the bodies in the inner binary 
has negligible mass compared to that of the outer perturber. Assuming, as done in previous treat- 
ments, that the z-component of the inner orbit's angular momentum (perpendicular to the invariable 
plane) is conserved can produce erroneous results for various astrophysical systems of interest. 
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1. INTRODUCTION 

Tri ple star systems are believed to bevery common 
(e.g., iTokovminl Il997t lEggleton et aTj 120071 ). From dy- 
namical stability arguments these must be hierarchi- 
cal triples, in which the (inner) binary is orbited by a 
third body on a much wider orbit. Pro bably more than 
50% of bright stars a re at least double (|TokoviniiJll997t 
lEggleton et al .112007ft . Given the selection effects against 
finding faint and distant companions we can be reason- 
ably confident that the propor tion is actually substan- 
tially greater. iTokovinml (|1997l ) showed that 40% of bi- 
nary stars with period < 10 d in which the primary is a 
dwarf (0.5—1.5 M@) have at least one additional compan- 
ion. He found that the fraction of triples and higher mul- 
tiples amo ng binaries with period (10 — 100 d) is ~ 10%. 
Moreover, iPribulla fe Rucinskil (12006ft have surveyed a 
sample of contact binaries, and noted that among 151 
contact binaries brighter than 10 mag., 42±5% are at 
least triple. 

Many close stellar binaries with two compact objects 
are likely produced through triple evolution. Secular ef- 
fects (i.e., coherent interactions on timescales long com- 
pared to the orbital per i od), and sp ecifically Kozai- 
Lidov cycling (|Kozailll962l lLidovlll962L see below), have 
been proposed as a n important element in the evolution 
of triple stars (e.g. lHarrington|[l 969: M azeh fc Shahaml 
1979t iKiseleva et all 119981 lFabrvckv fe Tremaind 120071: 
Perots fe Fabrv ckv 2009)7 In addition, Kozai-Lidov cy- 
cling has been suggested to play an important role in 
both the growth of black holes at the centers of dense star 
cluster s and the formation of short-period binary black 
holes rfWenl [2001 IMiller fc Hamilton! 120021: IBlaes et al.l 
12002ft . Recently, llvanova et all (|2010f ) showed that the 
most important formation mechanism for black hole 
XRBs in globular clusters may be triple-induced mass 
transfer in a black hole- white dwarf binary. 

Secular perturbations in triple systems also p lay an 
important role in planetary system dynamics. iKozail 
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(|1962ft studied the effects of Jupiter's gravitational per- 
turbation on an inclined asteroid in our own solar sys- 
tem. In the assumed hierarchic al configurat ion, treating 
the asteroid as a test particle, IKozail (|1962ft found that 
its inclination and eccentricity fluctuate on timescales 
much larger than its orbital period. Jupiter, assumed to 
be in a circular orbit, carries most of the angular mo- 
mentum of the system. Due to Jupiter's circular orbit 
and the negligable mass of the asteroid, the system's po- 
tential is axisymmctric and thus the component of the 
inner orbit's angular momentum along the total angu- 
lar mo mentum is conserved during the evolution. IKozail 
(|1979ft also showed the importance of secular interactions 
for the dynamics of comets (see alsolQiiinn et al.lfl99Ct 
iBailev et al1ll992HThomas fc Morbidellilll996f ). The evo- 
lution of the orbits of binary minor planets is dominated 
by the secular gravit ational perturbation from the sun 
(|Perets fc N aoz 2009J); properly accounting for the result- 
ing secular effects — including Kozai cycling — accurately 
reproduces the binary minor planet orbital distribu - 
tion seen tod ay dNaoz et al.l 120101 : j Grundv et all 12011ft. 
In 
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addi tion IKinoshita fc NakaiT (11991ft. . 
1999ft. ICarruba et all (j2002ft . iNesvornv et al.l (|2003ft . 

iCuk fc Burnl (f2004ft and IKinoshita fc Nakail ((2007ft sug- 
gested that secular interactions may explain the signifi- 
cant inclinations of gas giant satellites and Jovian irreg- 
ular satellites. 

Similar analyses ha ye been applied to the orbits of ex- 
trasq la r planets (e.g..llnnanen et al.lll997tlWu fc Murray 
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Maoz et al. | l201ll : lCorreia et al.ll201lft . iNaoz et all (|2011ft 
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considered the secular evolution of a triple system con- 
sisting of an inner binary containing a star and a Jupiter- 
like planet at several AU, orbited by a distant Jupiter-like 
planet or brown-dwarf companion. Perturbations from 
the outer body can drive Kozai-like cycles in the inner 
binary, which, when planet-star tidal effects are incorpo- 
rated, can lead to the capture of the inner planet onto 
a close, highly-inclined or even retrograde orbit, similar 
to the orbits of the observed retrograde "hot Jupiters." 
Many other studies of exoplanet dynamics have consid- 
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ered similar systems, but with a very distant stellar bi- 
nary companion acting as pcrturbcr. In such systems, 
the outer star completely dominates the orbital angular 
momentum, and the problem reduces to test-particle evo- 
lution in the lowest level of approximation, which leads to 
the conservation of the z -component of the inner orbit's 
angu l ar momentum (e.g. [Wu fe Murray! 120031: IWu et alj 
[2007t iFabrvckv fe Tremainell2007t iTakeda et al.ll2008t) . 

In early studies of high-inc lination secular perturba- 
tions (jKoza i 1962; Lidov 1961), the outer orbit was circu- 
lar and again dominated the orbital angular momentum 
of the system. In this situation, the component of the 
inner orbit's angular momentum along the z-axis is con- 
served. In many later studies the assumption that the 
z-component of the inner orbit's angula r momentum is 
const ant was built into the equations fe.g.lEggleton ct al.l 
19981 iMikkola, fc Tanikawal 119981 iZdziarski et al.1 120071 : 
Takeda et a l. 2008). In fact these studies are only valid 



in the limit of a test particle forced by a perturber on a 
circular orbit. To leading order in the ratio of semima- 
jor axes, the double averaged potential of the outer orbit 
is axisymmetric (even for an eccentric outer perturber), 
thus if taken to the test particle limit, this results in 
a conservation of the z-component of the inner orbit's 
angular momentum. We refer to this limit as the "stan- 
dard" treatment of Kozai oscillations, i.e. quadrupolc- 
level approximation in the test particle limit (test parti- 
cle quadrupole, hereafter TPQ). 

In this paper we show that a common mistake in the 
Hamiltonian treatment of these secular systems can lead 
to the erronious conclusion that the z-component of the 
inner orbit's angular momentum is constant outside the 
TPQ limit; in fact, the z-component of the inner orbit's 
angular momentum is only conserved by the evolution 
in the test-particle limit and to quadrupole order. To 
demonstrate the error we focus on the quadrupole (non- 
test-particle) approximation in the main body of the pa- 
per, but we include the full octupole-order equations of 
motion in an appendix. 

This paper is organized as follows. We first present the 
general framework (0; we then derive the complete for- 
malism for the quadrupolc-lcvcl approximation and the 
equations of motion (SJ3J), we also develop the octupole- 
level approximation equations of motion in §Q We dis- 
cuss a few of the most important implications of the cor- 
rect formalism in Sj5] We also compare our results with 
those of previous studies (S}6]) and offer some conclusions 

m m 

2. HAMILTONIAN PERTURBATION THEORY FOR 
HIERARCHICAL TRIPLE SYSTEMS 

Many gravitational triple systems are in a hierarchi- 
cal configuration — two objects orbit each other in a rel- 
atively tight inner binary while the third object is on a 
much wider orbit. If the third object is sufficiently dis- 
tant, an analytic, perturbative approach can be used to 
calculate the evolution o f the system. I n the usual sec- 
ular approximation (e.g.. iMarchaJ 119901 ) . the two orbits 
torque each other and exchange angular momentum, but 
not energy. Therefore the orbits can change shape and 
orientation (on timescales much longer than their orbital 
periods), but not scmimajor axes (SMA). 

We first define our basic notations. The system con- 
sists of a close binary (bodies of masses mi and 7712) and 
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Fig. 1. — Coordinate system used to describe the hierarchical 
triple system (not to scale). Here 'cm.' denotes the center of mass 
of the inner binary, containing objects of masses ra\ and 7712. The 
separation vector ri points from 777 1 to 7712; r2 points from 'cm.' 
to 7713. The angle between the vectors ri and r-2 is <&. 



a third body (mass 7773). It is co nvenient to describe 
the o rbits using Jacobi coordinates (jMurrav fc Dermottl 
120001 p. 441-443). Let ri be the relative position vector 
from m\ to mi and v 2 the position vector of 7773 relative 
to the center of mass of the inner binary (see fig. H}. Us- 
ing this coordinate system the dominant motion of the 
triple can be divided into two separate Keplerian orbits: 
the relative orbit of bodies 1 and 2, and the orbit of 
body 3 around the center of mass of bodies 1 and 2. The 
Hamiltonian for the system can be decomposed accord- 
ingly into two Keplerian Hamiltonians plus a coupling 
term that describes the (weak) interaction between the 
two orbits. Let the SMAs of the inner and outer orbits 
be a\ and a 2l respectively Then the coupling term in 
the complete Hamiltonian can be written as a power se- 
ries in the ratio of the semi-major axes a = ai/02 (e.g., 
lHarringtonlll968[ ) . In a hierarchical system, by definition, 
this parameter a is small. 

Th e complete Hamil tonian expanded in orders of a is 
(e.g.. lHarringtonlll968l L 
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where k 2 is the gravitational constant, Pj are Legendrc 
polynomials, $ is the angle between ri and r 2 (see Figure 
□J and 
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invariable plane 

Fig. 2. — Geometry of the angular momentum vectors. We show 
the total angular momentum vector (Gtot), the angular momentum 
vector of the inner orbit (Gi) with inclination i\ with respect to 
Gtot and the angular momentum vector of the outer orbit (G2) 
with inclination 12 with respect to Gtot- The angle between Gi 
and G2 defines the mutual inclination i tot = ii+12- The invariable 
plane is perpendicular to Gtot- 



Note t hat we have followed the convention of lHarringtonl 
(|1969t ) and chosen our Hamiltonian to be the negative of 
the total energy, so that H > for bound systems. 

We adopt the canonical variables known as Dclau- 
nay's elements, which provide a particularly convenient 
dynamical description of ou r three-body system (e.g. 
iValtonen fc Karttunenll2006D . The coordinates are cho- 
sen to be the mean anomalies, l\ and l 2 , the longitudes of 
ascending nodes, hi and h 2 , and the arguments of peri- 
astron, g\ and g 2 , where subscripts 1, 2 denote the inner 
and outer orbits, respectively. Their conjugate momenta 
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v //c 2 (mi + m 2 + m 3 )a 2 



Gi=LiJ\-e\ 1 G 2 = L 2 J1 



(3) 



(4) 



and 



Hi = G\ cosii , H 2 = G 2 cosi 2 



(5) 

Note that G\ and G 2 are also the magnitudes of the an- 
gular momentum vectors (Gi and G2), and Hi and H 2 
are the z-components of these vectors. Figure [2] shows 
the resulting configuration of theses vectors. The follow- 
ing geometric relations between the momenta follow from 
the law of cosines: 



cos i to t 
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H 2 = 



2Gtot 

r tot + ( - r 2 — U l 



2G t 



(6) 

(7) 
(8) 



where Gtot = Gi + G2 is the (conserved) total angu- 
lar momentum, and the angle between Gi and G2 de- 



fines the mutual inclination i to t = ii + *2- From eqs. ([7]) 
and ([8]) wc find that the inclinations ii and i 2 are deter- 
mined by the orbital angular momenta: 
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(10) 



In addition to these geometrical relations we also have 
that 

Hi+H 2 = G tot = const . (11) 

The canonical relations give the equations of motion: 
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(12) 
(13) 
(14) 



where j = 1,2. Note that these canonical relations have 
the opposite si gn relative to the usual relations (e.g., 
lGoldsteinlll950[ ) because of the sign convention we have 
chosen for our Hamiltonian. Finally we wri te the Hami l- 
tonian through second order in a as (e.g., lKozaillT962[ ) 
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where the mass parameters are 
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3. SECULAR EVOLUTION EQUATIONS TO QUADRUPOLE 

ORDER 

In this section, we derive the secular equations of mo- 
tion to the quadrupole-level, where in Appcndix[S]we de- 
velop the complete quadrupole-level secular approxima- 
tion. The main difference between the derivation shown 
here (see also Appendix |A| and those of pr evious stud- 
ies lies in the "elimina tion of nodes" (e.g., iKozail 119621 : 
Ueffervs fc Moserlll966D . This is related to the transition 
to a coordinate system with the total angular momentum 
along the z-axis, which is know n as the invariable plane 
(e.g.. iMurrav fc Dermottl 120001 ). In this coordinate sys- 
tem (see Figure [5]), the longitudes of the ascending nodes 
differ by it, i.e., hi — h 2 = tt. Conservation of angular 
momentum implies that this relation holds at all times. 
Many previous works have exploited it to explicitly sim- 
plify the Hamiltonian. However, as we explain below 
(and sec also %.!$ , this substitution leads to the incorrect 
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conclusion that Hi = H 2 = when the canonical equa- 
tions of motion are applied; thus, some previous stud- 
ies incorrectly concluded that the z-components of the 
orbital angular momenta are always constant. We will 
show that one can still use the (incorrect) Hamiltonian 
found in previous studies (e.g., IKozailfl962t IHarringtonl 
I1969T) as long as the evolution equations for the inclina- 
tions are derived from the total angular momentum con- 
servation, instead of using the canonical relations. Of 
course, the correct evolution equations can also be cal- 
culated from the correct Hamiltonian, which we derive 
in this section. 

We note that there are some other derivations of the 
secular evo lution equations that avoid the elimination of 
the nodes dFarago fc Laskarll2010HLaskar fc Boul l201d: 
IMardling||2010t IKatz fc Dong||2011i) . and thus do not suf- 
fer from this error. 

Previous studies made the substitution hi — hi = 7r 
directly in the Hamiltonian (see H6.ip . After the substi- 
tution, the Hamiltonian is independent of the longitudes 
of ascending nodes (h\ and h 2 ), and thus gives an evo- 
lution where both Hi and H2 are constant. The substi- 
tution hi — hi = 7r is incorrect at the Hamiltonian level 
because it unduly restricts variations in the trajectory of 
the system to those where 5hi = Sh2 (see Appendix IE)) . 
After deriving the equations of motion, however, we can 
exploit the relation hi — hi = Ah = n, which comes from 
the conservation of angular momentum and the fact that 
Gi + G2 = G to t = Gt tz. This considerably simplifies 
the equations. 

The secular Hamiltonian is given by the average over 
the rapidly-varying /1 and Z 2 in equation (|15p (Appendix 
E]for more details) 

U 2 = -/ { [1 + 3 cos(2i 2 )} ([2 + 3e 2 ] [1 + 3 cob(2*i)] (19) 
8 

+ 30e 2 cos(2 5 i) sin 2 (ii)) + 3 cos(2Aft.) [10e 2 cos(2fl>i) 
x {3 + cos(2ii)} + 4(2 + 3e 2 ) sin(ii) 2 ] sin 2 (i 2 ) 
+ 12{2 + 3e 2 - 5e 2 cos(2g x )} cos(A/i) sin(2ii) sin(2i 2 ) 
+ 120e 2 sin(ii) sin(2i 2 ) sin(2gi) s'm(Ah) 
- 120e 2 cos(ii) sin 2 (i 2 ) sin(2gi) sin(2A/i)} , 
where 
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(20) 



Making the usual (incorrect) transformation Ah — > n, we 
get the quadrupole-levcl Hamilto nian that has appe ared 
in many previous works (see, e.g. lFord et al.ll2~000bD : 

U 2 (Ah -)■ tt) = G 2 {(2 + 3e 2 ) (3cos 2 z tot - l) (21) 

+ 15e 2 sin 2 i t otcos(2gi)} , 

where we have set ii + i% = «tot- Because this Hamilto- 
nian is missing the longitudes of ascending nodes (hi and 
h 2 ), previous studies concluded that the z-componcnt 
(i.e., vertical) angular momenta of the inner and outer 
orbits (i.e., Hi and H 2 ) are constants. 

We use the canonical relations [equations (fT2")l ] in order 
to derive the equations of motion from the Hamiltonian. 
In our treatment, both Hi and H 2 evolve with time be- 
cause the Hamiltonian is not independent of hi and h 2 . 



From eq. (J7J), we see that 
Gi 



Hi = 



G t 



a G 2 ^ 

(-not 



(22) 



and from eq. (jlip we see that Hi = —H 2 . The 
quadrupole- level Hamiltonian does not depend on g 2 ; 
thus the magnitude of the outer orbit's angular momen- 
tum, G 2 , is constant, and therefore 



Hi = 



Ctot 



(23) 



From relations (J12H14J) we have Hi = dU/dhi, and Gi = 
dH/dgi. The former gives 

Hi = — 30C 2 e 2 sinz 2 sinj to t sin(2gii) . (24) 

and the latter evaluates to 

G, = -30G 2 e 2 sin 2 i tot sin(2 5l ) . (25) 

Employing the law of sines, Gtot/sini to t = Gi/sini 2 = 
G 2 /sinii, equation ((24)) can also be written as 



Hi = - 



Gi 

Gtnt 



30C 2 ei sin i to tsm(2gi) 



(26) 



which satisfies the relation in eq. (|2U)) . The evolution of 
the arguments of periapse are given by 

<7i = 6G 2 {^-[4cos 2 H ot + (5cos(2<?i) - 1) (27) 



x (1 - e{ - cos 2 itot)] 
and 



G 2 



-[2 + e 2 (3-5cos(2 ffl ))U , 



|^^[2 + e 2 (3-5cos(2 5l ))] 



(28) 



+ 7^[4 + 6e 2 + (5 cos 2 i tot - 3)(2 + e 2 [3 - 5cos(2 5 i)] 
G 2 

Previous quadrupole-level calculations that made the 
substitution error in the Hamiltonian lack the 1/G 2 term 
in the latter equation. The evolution of the longitudes of 
ascending nodes is given by 



in 



and 



3G 2 



Gi sinii 



h 2 = - 



3G 2 

G 2 sini 2 



{2 + 3e 2 - 5e 2 cos (2 5l )} sin (2i tot ) (29) 



{2 + 3e 2 -5e 2 cos(2 5 i)}sin(2i t ot). (30) 



Using the law of sines, Gi sinii = G 2 sini 2 , from which 
we get hi = h 2 , as required by the relation hi — h 2 = 7r. 
In many systems it is useful to calculate the time evolu- 
tion of the eccentricity, obtained through the following 
relation: 



dej dej &H 



(31) 



dt dGj dgi 

3 This conserved quantity is lost at higher orders of the approx- 
imation; see i|4]and Appendix [C] 
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In the quadrupole approximation e 2 = G 2 = (which is 
not the case at higher order in a; see Appendix [C]) . The 
eccentricity evolution for the inner orbit is given by 



Co 



1 



Gi 



-30ei sin 2 i tot sin(2#i) 



(32) 



Another useful parameter is the inclination, which can 
be found through the z-component of the angular mo- 
mentum: 



d(cosii) Hi Gi 



dt 



G\ G\ 



COS71 



(33) 



and similarly for ii (but note again that G 2 — to 
quadrupole order). In Appendix [B] wc show that the 
quadrupole approximation leads to well-defined mini- 
mum and maximum eccentricity and inclination. The 
eccentricity of the inner orbit and the inner (and mutual) 
inclination oscillate. We also demonstrate in Appendix ISI 
that our formalism gives critical initial mutual inclination 
angles for large oscillations of 39.2° < i to t < 140.8° in 
the test-particle limit and wit h nearly-zero initial inner 
eccentricity, in agreement with iKozal ([1962D . 

4. OCTUPOLE-LEVEL EVOLUTION 

In Appendix [C] we derive the secular evolution equa- 
tions to octupole order. Many previous octupole-order 
derivations provided correct secular evolution equations 
for at least some of the elements, in spite of using 
the elimination o f nodes substitu t ion a t the Hamil- 
tonia n level (e.g. iHarringtonl 1968 , 119691: iSidlichovskyl 



19831: iKrvmolowski fe Maze 



; 



199ft iFord et al.l 12000 



Blaes et all 12003 ILee fc PeaM2003bl: iThompsoplHSlC 

This is because the evolution equations for e 2 , g 2 , <7i and 
ei can be found correctly from a Hamiltonian that has 
had h\ and h 2 eliminated by the relation h\ — h 2 = it; 
the partial derivatives with respect to the other coordi- 
nates and momenta are not affected by the substitution. 
The correct evolution of Hi and H 2 can then be derived, 
not from the canonical relations, but from total angu- 
lar momentum conservation. We discuss in more details 
the comparison between this work and previous analyses 
in 3H1 The full consequences of allowing the Hi to be 
dynamical are explored here for the first time. 

The octupole-level terms in the Hamiltonian can be- 
come important when the eccentricity of the outer orbit 
is non-zero, and if a is not negligible. We quantifying 
this by considering the ratio between the octupole to 
quadrupole-level coefficients, which is 
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(34) 



where C3 is the octupole-level coefficient [cq. (|C1|) ] and 
C2 is the quadrupole-level coefficient [cq. ([2(3]) ]. Wc define 
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(35) 



which gives the relative significance of the octupole-level 
approximation. This parameter has three important 
parts; first the eccentricity of the outer orbit (e2), second, 
the mass difference of the inner binary (mi and 7712) and 



the SMA ratic0. In the test particle limit (i.e., mi 3> 7772) 
we find that cm is reduced to t he oc tupole coefficient 
introd uced in iLithwick fc Naozl (120111 ) and iKatz et all 
(I20T1 . 



f'2 



1 



(36) 



This coefficient measures the significance of the octupole 
contribution in the test particle limit. We will use here 
the general form (i.e., cm)- We label the behavior of a 
system for which ejvf <C 1 is not satisfied as "eccentric 
Kozai-Lidov" mechanism. 

The octupole terms vanish when e 2 = 0. Therefore if 
one artificially held e 2 = 0, in the test-particle limit the 
inner b ody's orbit w ould be given by the equations de- 
rived by|Kozai (1962]), i.e. by the test particle quadrupole 
equations. However, at octupole order the value of e 2 
evolves in time if the inner body is massive. Further- 
more, even if the inner body is masslcss, if the outer 
body has e 2 > then the inner body's behavior will 
also be different than in Ko zai's treatment. For e xam- 
ple. ILithwick fe Naozl (J2011D and IKatz etHU (J2011I ) find 
that the inner orbit can flip orientation (see below) even 
in the test-particle, octupole limit. The octupole-level 
effects can change qualitatively the evolution of a sys- 
tem. Compared to the quadrupole-level behavior, the 
eccentricity of the inner orbit can sometimes reach a 
much higher value. In some cases these excursions to 
very high eccentricities can be accompanied by a "flip" 
of the orbit with respect to the total angular momen- 
tum, i.e., starting with ii < 90° the inner orbit can 
eventually reach %i > 90° (see Figures [6] [9] for exam- 
ples) . Chjotk_bjhjv^r isalso_£C^shblj_^at_Jhe octupole 
level (|Lithwick fe Naozll201ll : IKatz et al.ll2011l) . but not 
at the quadrupole-level (see Appendix |BJ) . In contrast 
to the octupole-level behavior, the quadrupole-level ap- 
proximation leads to regular cycles in eccentricity and 
inclination, with well-defined maximum and minimum 
values, and it cannot produce flips for the inner orbit 
(again, see Appendix IB)) . 

Given the large, qualitative changes in behavior mov- 
ing from quadrupole to octupole order in the Hamil- 
tonian, is it possible that similar changes in the sec- 
ular evolution may occur at even higher orders? In- 
tuitively, the answer to this question lies in the elimi- 
nation of G 2 as an integral of motion at octupole or- 
der, leaving only four integrals of motion: the energy of 
the system, and the three components of the total an- 
gular momentum. There are no more integrals of mo- 
tion to be eliminated, and thus one might expect no 
more dramatic changes in the evolution when moving 
to even higher orders. It is possible to see this quantita- 
tively for specific initial conditions through comparisons 
with direct 77-body integrations. We compare our oc- 
tupole equations with direct 77- body integrations, using 
the M ercury software package (jChambers fe Migliorinil 
11997ft . We used both Bulirsch-Sto er and symplectic inte- 
grators (jWisdom fc Holmanlll991l ) and found consistent 
results between the two. We present the results of a 
typical integration compared to the integration of the 

4 Note here that the subscripts "1" and "2" refer to the inner 
bodies in mi and m^, but the subscript "2" refers to the outer 
body in e2- 



Naoz ct al. 



180 




2.5x10' 



Fig. 3. — Comparison between a direct integration (using a B- 
S integrator) and the octupolc-level approximation (see Appendix 
let . The red lines are from the integration of the octupolc-level 
perturbation equations, while the blue lines are from the direct 
numerical integration of the three-body system. Here the inner bi- 
nary contains a star of mass 1 Mq and a planet of mass 1 Mj , while 
the outer object is a brown dwarf of mass 40 Mj. The inner orbit 
has a\ = 6 AU and the outer orbit has a,2 = 100 AU. The initial 
eccentricities are e\ = 0.001 and e^ = 0.6 and the initial relative 
inclination itot = 65°. The thin horizontal line in the top panel 
marks the 90° boundary, separating prograde and retrograde or- 
bits. The initial mutual inclination of 65° corresponds to an inner 
and outer inclination with respect to the total angular momentum 
(parallel to z) of 64.7° and 0.3°, respectively. The arguments of 
periccnter of the inner and outer orbits are initially set to zero. The 
SMA of the two orbits (not shown) are nearly constant during the 
direct integration, varying by less than 0.02 percent. The agree- 
ment in both period and amplitude of oscillation between the direct 
integration and the octupolc-level approximation is quite good. 



octupolc-level secular equations in Figure [3] The ini- 
tial conditi o ns (se e caption) for this system are those of 
iNaoz et ail ()2011h . Figure 1. We find good agreement 
between the direct integration and the secular evolution 
at octupole order. Both show a beat-like pattern of ec- 
centricity oscillations, suggesting an interference between 
the quadrupole and octupole terms, and both methods 
show similar flips of the inner orbit. 

5. IMPLICATIONS 

The iKozail (|1962h and iLidovl (|1962t l equations of mo- 
tions are correct to quadrupole order and for a test par- 
ticle. Using them outside this limit can lead to incorrect 
results. Here we discuss a few examples and the im- 
portance of higher-order effects in various astrophysical 
settings. 

5.1. Massive Inner Object at the Quadrupole Level 

The danger with working in the wrong limit is apparent 
if we consider an inner object that is more massive then 
the outer object. While the standard formalism incor- 
rectly assumes that the orbit of the outer body is fixed 
in the invariable plane, and therefore the inner body's 
vertical angular momentum is constant, the quadrupolc- 
level equations presented in Section [3] do not. 

We compare the two formalisms in Figure |4j We con- 
sider the triple system PSR B1620— 26 located near the 



core of the globular cluster M4. The inner binary con- 
tains a millisecond radio puls ar of mi = lAMp , and a 
comp anion of m 2 = 0.3 M^ (jMcKenna fc Lvnel [l98cl . 
From iFord et aLl (j2000af ) , we adopt parameters for the 
outer perturbcr of m^ = 0.01 M Q and e 2 = 0.45. The 
inner binary separation is a\ — 5 AU while a 2 = 50 AU. 
We initialize the system with i tot = 70°, g\ = 120° and 
g-2 = 0° and ei = 0.5. Note that the actual measured in- 
ner binary eccentricity is e\ ~ 0.045, however in order to 
illustrate the difference we adopt a higher value. There is 
no logical reason to assume the observed eccentricity as 
the initial eccentricity when modeling the formation of a 
system, since different physical processes can contri bute 
to ec centricity damping (for example, tidal fr iction (iHutl 
Il980f ) and mass transfer (|Sepinskv et al.ll20 10)). The ini- 
tial mutual inclination of 70° corresponds to an inner and 
outer inclination with respect to the total angular mo- 
mentum (parallel to z) of 6.75° and 63.25°, respectively. 
Remember that we consider the evolution of the system 
to quadrupole order for comparison, even though there is 
no apriori reason to truncate the evolution at this order, 
especially since cm = 0.036. We have verified, however, 
that octupole order effects do not qualitatively change 
the behavior. This is because the outer companion mass 
is low, and hence the inner orbit does not exhibit large 
amplitude oscillation^- 

For the comparison, we do not compare the (constant) 
Hi from the TPQ formalism to the (varying) Hi of the 
correct formalism. Instead, we compare the (constant) 
Hi from the TPQ formalism with Gi cos i to t, which is the 
vertical angular momentum that would be inferred in our 
formalism if the outer orbit were instantaneously in the 
invariable plane, as assumed in the TPQ formalism. 

In Figure [4j the mutual inclination oscillates between 
106.7° to 57.5°, and thus crosses 90°. These oscillations 
arc mostly due to the oscillations of the outer orbit's 
inclination, while ii does not change by more than ~ 1° 
in each cycle. Clearly, the outer orbit does not lie in the 
fixed invariable plane! Figure [4j bottom panel, shows 

yl — ef cos itot) which, in the TPQ limit, is the vertical 
angular momentum of the inner body. 

We can evaluate analytically the error introduced by 
the application of the TPQ formalism to this situation. 
We compare the vertical angular momentum (Hi ) as cal- 
culated here to H 1 y = L\\J\ — e\ cosi = const.. The 

relative error between the formalisms is H x /Hi — 1. 
In Figure [5] we show the ratio between the inner or- 
bit's vertical angular momentum in the TPQ limit (i.e., 
H ± ^ = Gi cosi) and equation (|2"4"|) as a function of the 
total angular momentum ratio, G1/G2, for various incli- 
nations. Note that this error can be calculated without 
evolving the system by using angular momentum conser- 
vation, equation ([6]). The TPQ limit is only valid when 
G1/G2 < 10- 4 . 

5.2. Planetary Dynamics 

Recent measurements of the sky-projected angle be- 
tween the orbits of several hot Jupitcrs and the spins 

5 Unlike the test particle octupol e-level approximation 
IILithwick fc Naozl [2011] ; IKatz et al.l [20liT ). backreaction of the 
outer orbit may suppress the eccentric Kozai effect. We address 
this in further detail in Teyssandier et al. in Prep. 
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Fig. 4. — Comparison between the standard TPQ formalism 
(dashed blue lines) evolution and our methods (solid red lines) 
for the case of PSR B1620— 26. Here the inner binary is a mil- 
lisecond pulsar of mass 1.4 Mq with a companion of m,2 = 0.3Mq, 
and the outer body has mass 7713 = O.OIMq. The inner orbit ha s 
ai = 5 AU and the outer orbit has 02 = 50 AU I IFord et alJ2000af l. 
The initial eccentricities are e\ = 0.5 and e2 = 0.45 and the ini- 
tial relative inclination itot = 70°. The thin horizontal line in the 
top panel marks the 90° boundary, separating prograde and ret- 
rograde orbits. The initial mutual inclination of 70° corresponds 
to an inner and outer inclination with respect to the total angu- 
lar momentum (parallel to z) of 6.75° and 63.25°, respectively. 
The argument of pericenter of the inner orbit is initially set 120°, 
while the outer orbit's is set to zero. We consider, from top to 
bottom, the mutual inclination itot > the inner orbit's eccentricity 

and Jl-ej cos itot, which the standard formalism assumes to be 

constant (dashed line). 
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Fig. 5. — The ratio between the correct, changing z-component 
angular momentum, Hi, and the assumption often used in the liter- 
ature which is H 1 y = Gi cos i. This ratio was calculated analyti- 
cally for various total angular momentum ratios, G1/G2, and incli- 
nations. The curves, from bottom to top, have i = 40, 50, 60, 70, 80 
and 89 degrees. 



Fig. 6. — Evolution of a planetary system with m\ = IMq, 
m,2 = lMj and 7713 = 2Mj. with ai = 4 AU and 02 = 45 AU. 
We initialize the system at t = with ei = 0.01, e2 = 0.6, g\ = 
180°, <?2 = 0° and itot = 67°. For these initial conditions ii = 
57.92° and i2 = 9.08°. The z-componcnts of the orbital angular 
momenta, H\ and H2, arc shown normalized to the total angular 
momentum of each orbit. The inner orbit flips quasi-pcriodically 
between prograde (ii < 90°) and retrograde (ii > 90°). 



of their host stars have shown that roughly one in 
four is retrograde (jGaudi fc Winnl l2007t iTriaud et al.l 
2010). If these planets migrated in from much larger 
distances through their interaction with the protoplane- 
tary disk (jLin fc Papalo izou 1986; Massc t fc Papaloizoul 
120031 ). their orbits should have low eccentricities and 
mclinationio Disk migration scenarios therefore have 
difficulty accounting for the observed retrograde hot 
Jupiter orbits. An alternative migration scenario that 
can account for the retrograde orbits is the secular in- 
teractio n between a planet and a binary stellar com- 
panion (IWu fc Murravl2003HFabrvckv fc Tremamell2007l ; 
IWu et al.H2007HTakeda et al.ll2008t ICorreia et al.H20l¥ . 
For a very distant companion (em < 1) the quadrupole 

test-particle approximation applies, and yl — ef cosii 
is nearly constant. Although this forbids orbits that are 
truly retrograde (with respect to the total angular mo- 
mentum of the system), if the inner orbit begins highly 
inclined relative to the outer star's orbit and aligned with 
the spin of the inner star, then the star-planet spin-orbit 
angle can change by more than 90° during the secular 
evoluti on of the system, producing apparently retrograd e 
orbits (jFabrvckv fe Tremaindl20Cm ICorreia et al.ll2011l) . 
Nonetheless, a difficulty with this "stellar Kozai" mech- 
anism is that even with the most optimistic as sumptions 
it ca n only produce < 10% of hot Jupiters (IWu et al.l 

120071). 

iNaoz et"ail ()201lD considered planet-planet secular in- 
teractions as a possible source of retrograde hot Jupiters. 
In this situation em is not small, requiring computation 

6 This assumption can be invalid if there are significant mag- 
netic interactio ns between the star and the protoplanetary disk 
IILai et al.H2010l) or if there is an episode of planet-planet scattering 
follow ing plan et formation ( Chattcrjcc et al. 2008; Nagasa wa et ail 
I2008T) see also lMerritt et all J20091) . 
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Fig. 7. — Zoom-in on part of the evolution of the point-mass 
planetary system in Figure[6] In this zoom-in, we can see that flips 
in the inner orbit — i\ crossing 90° — are associated with excursions 
to very high eccentricity. 



of the octupole- level secular dynamics. In Figures[B]and[7] 
we show the evolution of a representative configuration, 
where mi = 1 M Q , m 2 = lMj and 7713 = 2Mj, with 
Oi = 4 AU and a 2 = 45 AU. For this configuration, 
cm = 0.083. Flips of the inner orbit are associated with 
evolution to very high eccentricity (see Figure [7]). 

5.3. Solar System Dynamics 

iKozail (119621 ) studied the dynamical evolution of an 
asteroid due to Jupiter's secular perturbations. He as- 
sumed that Jupiter's eccentricity is strictly zero. How- 
ever, Jupiter's eccentricity is ~ 0.05, and thus study- 
ing the evolution of a test particle in the asteroid belt 
(ai ~ 2 — 3 AU) places the evolution in a regime where 
the eccentric K ozai-Lidov effect could be significant, wit h 
e M = e = 0.03 (jLithwick fe Naoz|[20TTI : lKatz et al.ll20TTI ). 

We considered the evolution of asteroid at 2 AU (as- 
sumed to be a test particle) due to Jupiter at 5 AU 
with eccentricity of e 2 = 0.05. We set ei = 0.2, 
*tot = 65° and g\ = g 2 = 0° initially. The aster- 
oid is a test particle and therefore i\ ss itot- I n Fig- 
ure [5] we compa re the evolution of an asteroid using the 
TPQ limit fe.g .IKozailll962t lTho mas fc Morbideli1 [l99l 
iKinoshita fc NakaH 120071 ) and the octupole-level evolu- 
tion discussed here. For this value of e, the eccentric 
Kozai-Lidov effect significantly alters the evolution of the 
asteroid, even driving it to such high inclination that the 
orbit becomes retrograde. Though we deal only with 
point masses in this work, note that the eccentricity is 
so high that the inner orbit's pericenter lies well within 
the sun. 

The value of e here is mainly due to the relative high 
a in the problem. The system is very packed which raise 
questions with regards to the actual validity of the ap- 
proximation in that regime. In fact, such high eccentrici- 
ties drive the pericenter of the asteroid to collide with the 
sun and the apo-center of the asteroid to approach about 
1 AU from Jupiter's orbit. To address this question we 
run a iV-body simulation using Mercury software package 



Fig. 8. — Evolution of an asteroid due to Jupiter secular grav- 
itational perturbations (Kozai 1962). We consider mi = IMq, 
m-2 — > and 1713 = 1 Mj, with a\=2 AU and 02 = 5 AU. We ini- 
tialize the system at t = with e\ = 0.2, e2 = 0.05, 91 = 92 = 0° 
and «tot = 65°. We show the TPQ limit evolution (cyan lines) 
and the octupole equations (red lines). The thin horizontal dotted 
line in the top panel marks the 90° boundary, separating progradc 
and retrograde orbits. The inner orbit flips quasi-periodically be- 
tween prograde (i\ < 90°) and retrograde (i± > 90°). We also 
show the result of an JV-body simulation (blue lines). The thin 
horizontal dotted line in the bottom panel marks the solar radius 
as 1 — ei = Rq/cii. 



(|Chambers fc Migliomul 119971 ). We used both Bulirsch- 
Stoer and symplectic integrators (|Wisdom fc Holmanl 
119911 ). The results are depicted at Figure which show 
that the TPQ limit is indeed in adequate for the sys- 
tem. In addition the octupole-level approximation has 
some deviations from the direct TV-body integration, and 
does not follow the direct integration results in the high 
eccentricity regime. Note that the evolution of the as- 
teroid in the direct integration had resulted with a col- 
lision with the Sun. In reality, highly eccentric asteroids 
do not live very long in the solar system, also due to 
planet-crossing. We also note that the assuming zero 
eccentricity for Jupiter results in consistent results (we 
tested the case for a\ — 2 AU) be tween the secular evo- 
lution and the direct integration (jThomas fc Morbidellil 
119961 see also). Note also that Kozai mentions that the 
TPQ limit may not be correct, both due to the relatively 
large value of a and the non-negligable eccentricity of 
Jupiter, but proceeds anyway with the theory because it 
is analytically tractable. 

5.4. Triple Stars 

The evolution of triple stars has been studied by 
many auth ors using the standard ( TPQ) formal- 



ism 



1998: 
1998' 



1 



iMazeh fc Shahaml 119791: lEggleton et al 



Kiseleva et al 



1998 



Egglcton & Kisclcva-Egglcton 



Fabrvckv fc Tremaind I2007t IPerets fc Fabrvckvl 




In some cases the corrected formalism derived here can 
give rise to qualitatively different results. We show that 
some of the previous studies should be repeated in order 
to account for the correct dynamical evolution, and give 
one example where the eccentric Kozi-Lidov mechanism 
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Fig. 9. — An example of dramatically different evolution between 
the quadruple and octupole approximations for a triple-star sys- 
tem. The system has mi = IMq, rri2 = O.IMq and 7713 = OAMq, 
with a\ = 2 AU and 02 = 11 AU. We initialize the system with 
ei = 0.01, e 2 =0-6, .91 = 1 4 5°, g 2 = 0° and i t ot = 65° accord- 



ing to Fabrycky ik, Trcmainc (2007) Monte-Carlo simulations. For 
these initial conditions i\ = 58.1° and 12 = 6.9°. We show both 
the (correct) quadrupolc-lcvel evolution (light-blue lines) and the 
octupole-level evolution (red lines). H± and H2, the 2-components 
of the angular momenta of the orbits, are normalized to the total 
angular momentum. Note that the octupole-level evolution pro- 
duces periodic transitions from prograde to retrograde inner orbits 
(relative to the total angular momentum), while at the quadrupolc- 
level the inner orbit remains prograde. 



d ramatically changes the evolu tion. 

iFabrvckv k Trema inc (2007|) studied the distribution 
of triple star properties using Monte-Carlo simulations. 
We choose a particular system from their triple-star suite 
of simulations to illustrate how the dynamics including 
the octupole order can be qualitatively different from 
what would be seen at quadrupole order, when em is 
not negligible^. We adopt the following parameters: 
mi = 1M Q , m 2 = 0.1M Q and m 3 = 0.4M Q , ai = 2AU 
and a 2 = 11 AU. We initialize the system at t — with 
d = 0.01, e 2 = 0.6, g x = 145°, g 2 = 0° and z tot = 65°, 
corresponding to i\ = 58.1°, z 2 = 6.9° and cm = 0.14. 
The evolution of the system is shown in Figure [9] At 
octupole order, the inclination of the inner orbit oscil- 
lates between about 40° and 140°, often becoming ret- 
rograde (relative to the total angular momentum) , while 
the quadrupole-order behavior is very different and the 
inner orbit remains always prograde. The octupole-ordcr 
treatment also gives rise to much higher eccent ricities 
(jKrvmolowski fc Mazebl [l999t iFord etaLI l2000b[) . The 
evolution shown in Figure |H] is for point-mass stars; in 
reality, these high-eccentricity excursions would actually 
drive the inner binary to its Roche limit, leading to mass 
transfer. 

The possibility of forming blue stragglers through secu- 
lar interactions in tripl e star s ystem s has been suggeste d 
by iPerets k Fabrvckvl (120091) and iGeller et al.l l20ll 



As shown in iKrvmolowski k Mazehl ( 19991 ); IFord et all 



(|2000bl) and in the example above the minimum peri- 
center that the inner binary can reach can vary (from 
few percents to orders of magnitude depending on em). 
Thus, it suggests that the correct formalism may increase 
the likelihood of such a formation mechanism for blue 
stragglers. 

For many years CH Cygni was considered to be an 
interesting triple candidate becau se it exhibits two clear 
distin gui shable perio ds (e.g . iDonnison k Mikulskis 
ISkopal et al.l J1998T; IMikkola k Tanikawa 
However, a triple sys- 
the TPQ Kozai mechanism 



[1993 ) 



on 



19981: iHinkle et a" 
tern model based 

(|Mikkola k Tanikawalll998T ) did not reproduce t he ob - 
served masses of the system (|Hinkle et al.l 119931 I2009D . 
Applying the corrected formalism in this paper to the 
system parameters derived in IMikkola k Tanikawa! 
(|1998|) gives a very different evolution than in the TPQ 
formalism^. Therefore, it seems likely that an analysis 
based on the formalism discussed in this paper would 
give a significantly different fit. In Figure[TU]wc illustrate 
the differences between the TPQ, correct quadrupole, 
and octupole evolution of the s ystem. The best-fit 
param eters of the system from IMikkola k Tanikawal 
(|1998[) are as follows: mi = 3.51M , m 2 = 0.5M Q and 
m 3 = O.9O9M , oi = 0.05 AU and a 2 = 0.21 AU. We 
initialize the system at t = with e\ = 0.32, e 2 = 0.6, 
.gi = 145°, g 2 = 0° and itot = 72°, corresponding to 
z'l = 57.01°, i 2 = 14.98° and e m = 0.14. We allowed for 
a freedom in our choice of e 2 , g\ , g 2 and itot since the the 
best fit was found using the TPQ limit, at which e 2 is 
fixed. Note that the choice of the inner eccentricity does 
not strongly influence the evolution while the choice of 
the outer orbit's eccentricity does. 

6. COMPARISON WITH PREVIOUS STUDIES 



Kozai (1962) studied the motion of an inclined aster- 
oid due to perturbations from Jupiter. He derived the 
Hamiltonian for this system to high order in a , assum- 
ing a circular orbit for Jupiter. He then truncated the 
expansion at quadrupole order in a to derive the secular 
evolution equations for the asteroid; his equations thus 
correctly describe the test-particle quadrupole (TPQ) 
limit. However, Kozai's equations were later applied in- 
correctly in other studies. Kozai's equations imply that 
Hi = const, but outside the test-particle limit in the 
quadrupole-order evolution we have seen that Hi is no 
longer constant. Moreover, even when em is small, the 
octupole-order effects can lead to qualitatively different 
orbits than predicted at quadrupole order. 

6.1. Elimination of the Nodes and the Problem in 
Previous Quadrupole- Level Treatments 

Since the total angular momentum is conserved, the 
ascending nodes relative to the invariable plane follow a 
simple relation, hiit) = /i 2 (£) — tt (see Appendix iDl) . If 
one inserts this relation into the Hamiltonian, which only 
depends on hi — /i 2 , the resulting "simplified" Hamilto- 
nian is independent of hi and /i 2 . One might be tempted 
to conclude that the conjugate momenta Hi and H 2 are 



' The extent of the Fabrvckv & Trcmainc (2007) phase space 
over which ejvf is not negligible requires further investigation. 



8 Mikkola & Tanikawa ( 1998) also found somewhat different set 
of parameter when producing a fit for data set with less weight 
for the data of 1983 due to large noise in the active phase of the 
system. 
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Fig. 10. — An example of dramatically different evolution 
between the quadruple and octupole approximations for a 
triple star system represen ting the best-fit parameters from the 
IMikkola fc Tanikawal 1 119981) analysis of CH Cygni. The system 
has m 1 = 3.51M , m 2 = O.5M and m 3 = 0.909M Q , with 
a\ = 0.05 AU and a.2 = 0.21 AU. We initialize the system with 
ei = 0.32, e 2 = 0.6, g\ = 145°, gi = 0° and t to t = 72°. For these 
initial conditions ii = 57.02° and %i = 14.98°. We show both the 
(non-TPQ) quadrupole-level evolution (light-blue lines) and the 
octupole-levcl evolution (red lines). Hi and H2, the z-components 
of the angular momenta of the orbits, are normalized to the total 
angular momentum. Note that the octupole-level evolution pro- 
duces periodic transitions from prograde to retrograde inner orbits 
(relative to the total angular momentum), while at the quadrupole- 
level the inner orbit remains prograde. To avoid clatter in the figure 
we omitted the TPQ limit result, we note however, that the evo- 
lution of the inclination and eccentricity are similar to the general 
quadrupole-level approximation, but with constant H±2- 



constants of the motion. However, that conclusion is 
false. This incorrect argument has been made by a num- 
ber of authorfl 

We now show that H\ and H 2 are constant only in 
the TPQ limit. Outside of that limit they are not con- 
stant, and taking them to be constant breaks the con- 
servation of the total angular momentum of the system. 
The quadrupole-level Hamiltonian (Eq. [T9|) depends on 
the angles as follows: H = H(gi,hi — h 2 ); that is, it 
is independent of 172 and only depends on the difference 
between h\ and h 2 . Using the incorrect elimination of 
the nodes discussed above, one would conclude that Hi 
is constant. Also, because the Hamiltonian is indepen- 
dent of #2, G2 = const. From the geometric relation in 
equation ©, Hi = (G t 2 ot + G? - G 2 2 )/(2G tot ), and the 
constantcy of the total angular momentum, Gtot , we have 
that G\ = const, but this is inconsistent with the depen- 
dence of T-L on 31. The error comes from assuming that 
H\ = const; in fact, at the quadrupole level, we have 



Hi - — — Gi, 

(-"•tot 



(37) 



which is consistent with both the geometric relation in 



9 For example, [Kozai ( 1962, p. 592) incorrectly argues that "As 
the Hamiltonian F depends on h and h' as a combination h — h' , 
the variables h and h' can be eliminated from F by the relation 
(5). Therefore, H and H' are constant." 



equation (J7J and the dynamical equations (|24[) and 
When Gi/Gtot -> 0, or, equivalently G 2 > Gi, H 1 -> 0, 
and the TPQ result is achieved 1 !. But, this is precisely 
the TPQ limit, where the outer orbit dominates the an- 
gular momentum of the system and therefore lies in the 
invariable plane. 

In general, using dynamical information about the 
system — in this case that angular momentum is con- 
served, implying that Gi + G2 = Gtot at all times and 
therefore hi — h 2 = it — to simplify the Hamiltonian is not 
correct. The derivation of Hamilton's equations relies 
on the possibility of making arbitrary variations of the 
system's trajectory, and such simplifications restrict the 
allowed variations to those which respect the dynamical 
constraints. Once Hamilton's equations are employed to 
derive equations of motion for the system, however, dy- 
namical information can be employed to simplify these 
equations. See Appendix [El for further discussion of this 
point. 

In our particular case, quations of motion for compo- 
nents of the system that do not involve partial deriva- 
tives with respect to hi or h 2 will not be affected by 
the node-elimination substitution. For this reason, it is 
correct to derive equations of motion for all components 
except for Hi and H 2 from the node-eliminated Hamil- 
tonian; expressions for Hi and H 2 can then be derived 
from conservation of angular momentum. This approach 
has been employed in at least one computer code for 
octupole evolution, though the discussion in the corre- 
sponding pap er incorrectly elim inates the nodes in the 
Hamiltonian (jFord et a l. 2000b). 

In some later studies . (ISidlichovskvl 

Innanen et al.l Il997t iKiseleva et al. 

Eggleton et all 119981: IMikkola k Tanikawal 



1983 



1998 



1998 



Kinoshita k Nakai 119991: lEggleton fc Kiseleva-Eggleton 



Wu k Murray! 120031: IValtonen fc Karttunen 
200a IFabTvckyfc Tremaind 120071: IWu et all 120071: 
Zdziarski et all 120071: iPerets k Fabrv ckv 200^7 the as- 



sumption that Hi = const was built into the calculations 
of quadrupole-level secular evolution for various astro- 
physical systems, even when the condition G2 3> Gi 
was not satisfied. Moreover many p r evious studies 
(ISidlichovskvl 119831: llnnanen et all 119971: IKiseleva et all 
19981: lEggleton et all 119981: IMikkola k Tanikawal 1199 



Kinoshita k Nakall 119991 : lEggleton k Kiseleva-Eggleton 



2001 1 : IWu fc Murravl 12003 : IValtonen k Karttunen 
mia iFabrvckv k Tremamel 120071 : IWu et all 120071 : 
Perets fc Fabrvckv! 120091 ) simply set «2 = 0. In fact, 
given the mutual inclination i, the inner and outer 
inclinations ii and ii are set by the conservation of total 
angular momentum [see equations (J9j) and (|10[) ]. 

In Figurc[5]we show the ratio between the inner orbit's 
vertical angular momentum in the TPQ limit (H 1 = 
Gicosi) and the Hi in the correct quadruple-level ap- 
proximation from the derivation shown here as a func- 
tion of the ratio G1/G2. From this figure, it is clear that 
the standard formalism is only valid in the TPQ limit, 
where G1/G2 < 10~ 4 (depending slightly on the mutual 

10 This is true as long as Gi is not larger then Gi/Gtot- How- 
ever, Gi is proportional to C2 which is in turn proportional to 
1/G§ ~ 1/G? 01 
limit. 



Therefore, GiGi/Gtot — > in the test- particle 
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180 




2.5x10' 



Fig. 11. — Comparison between direct integration and the 
octupole-level approximation. We consider both the wrong sign 
(i.e., G's > for mi > mi) and the correct sign (i.e., C3 < for 
mi > mi) for the octupole term. The dotted curve represents 
the result from the integration with C3 > 0, while the thick solid 
curve depicts the integration with the C3 < 0. The thin curve, 
overlapping the thick curve shows the result of direct numerical in- 
tegration of the three-body system using Burlish-Stoer integrator. 
We have used here the initial conditions in Figure [3] 

inclination). 

The quadrupole-level equations appearing com- 
monly in the literature are correct only in the 
limit of qu adrupole-level test particle (jKozail Il962t 
iLidovl Il962f ). Nevertheless, these equations have 
been applied outside this limit, to massive bodies in 

the inner binary i n many studi es in the literature 

(ISidlichovskvl 119831: llnnanen et all Il997t iKiseleva et all 



19981: lEggleton et al.l 7 l998t iMikkola fc Tanikawal 1199? 
Kinoshita fc Nakal 119991 : lEggleton fc Kiseleva-Eggleton) 
2001 1 : IWu fe Murravl 12003 iValtonen fe Karttunenl 
200l Ifabrvckv k Tremaind 120071: IWu et all l2007t 
Zdziarski et al.ll2007HPerets fc Fabrvcky||2009D 



6.2. Octupole-level Approximation and Truncation to 
Quadrupole 



The octupole-level Hamiltonian 

of mot ion we re previously derived 

" 19681 [19691): ISidhchovskvi (fl98l: IMarchall 



Krvmolowski fc Mazehl (119991): IFord et all 



and equations 

by lHarringtor 

199C 



Blaes et all (|2002H and ILee fc Peald (j2003bl) . Most of 



the equations of motion can be derived correctly when 
applying the elimination of the nodes — only the Hi and 
H2 equations are affected. These authors calculated 
the time evolution of the inclinations (i.e. H\ and H2 
from the total (conserved) angular momentum, and 
thus avoided the problem that arises when eliminating 
the nodes from the Hamiltonian. In appendix [C] we 
show the complete set of equations of motion for the 
octupole-level approximation, derived from a correct 
Hamiltonian, including the n odal terms. 

As prev iously note d by iBlaes et al.l (|2002l equa- 
tion 24), IFord et ail (|2000bl ) introduced a sign er- 
ror in the o ctupole coefficient , C3, which was later 
corrected in IFord et ail (|2004f ). The same sign er- 



ror also exists in Marchal] (1199 
cq. 17), iKrvmolowski fc Mazehl 



lLaskar fc Bouel (120101 their equation for T . 
ilfjh and ILee ~ 



Oh. ISidlichovskvl (fl983l 
(|l999l ea. t" 

(0,0) 



and 



eq. 6b) 

;. W e note 
that iThompsonl (|2010l ) and ILee fc Peald (|2003bl |ah used 
the correct sign. To settle this point we use the direct 
N-body simulation shown in Figure [3] and compare it 
to the integration of the octupole-level approximation 
equations (see Appendix [C]) with and without the minus 
sign. We show our r esul ts in Fig ure ITT1 and f i nd th at C3 
as defined here [eq. (|C1|) ] and in lBlaes et al.l <|2002f ) is in 
agreement with the direct N-body, and thus indeed this 
is the correct sign. We note that this sign error can be 
easily resolved if 52 —> .92 + tt (which may imply to the 
source of this confusion). 

As displayed here the octupole-level approximation 
gives rise to a qualitatively different evolutionary 
behavior for cases where €m [see eq. ([33])] is not negli- 
gible. We note that many previous studies applied the 
quadrupole-level approximation , which may lead to sig- 
nificantly different r esults (e.g.. iMazeh fc Shahamlll979 
QuinnetaTI 119901: iBailev et all 119921: llnnanen et alj 



19971: lEgg eton et all Il998t IMikkola fc Tanikaw; 



1998; 



Valtonen fc Karttunenl 



Eggleton fc Kiseleva-Eggleton! 



2007F 



Wuetal 



1J I2007t 

312Q09D 



2001; 



2006T iFabrvckv fc Tremaind 
iZdziarski etafl 120071: 
Perets fc Fabrvckvl 120091 ) . Neglecting the octupole- 
level approximation can cause changes in the dynamics 
varying from a few percent to completely different 
qualitative behavior. 

Some other derivations of octupole-order equations 
of motion dealt with the secular dynamics in a gen- 
eral way, without using Hamiltonian perturbation the- 
ory or elimination of the nodes dFarago fc Laskar 



201(1 lLaskar fc Boudl2Mol: iMardling 20101: IKatz fc"Poiig 



20111 ). In these works there were no references to the 
discrepancy between these derivations an d the previous 
studie s. Also, note that the results of iHolman et all 
(|1997l) are based on a direct N-body integration, and 
thus are not subject to the errors mentioned above. 

6.3. Comparisons with Specific Papers 

Many previous studies applied the test particle 
quadrupole-level equations (TPQ) in various astrophysi- 
cal settings, even in situations where those equations are 
not strictly applicable. We address some of these studies 
here in more detail. ^_^^ ^_^ 

As discussed above, iKozail (|1962l ) studied the effect of 
Jupiter perturbations to an inclined asteroid. He specif- 
ically assumed that Jupiter's eccentricity is zero. How- 
ever, as shown in Figure 8, taking into account Jupitcrs 
eccentricity (~ 0.05), produces a dramatically different 
evolution ary behavior, including retrog rade orbits for the 
asteroid. iThomas fc Morbidellil (| 19961 ) applied the same 
limit to the asteroid- Jupiter setting (see for example their 
Figure 2 for a 1 = 3 AU, where they explicitly show 
the (wrong) vertical angu lar momentum conservation). 
iKinoshita fc Nakai (J2007D deve loped an analytical solu- 
tion f or the TPQ limit (see also IKinoshita fc NakaillT99lL 
I1999D . however, they have applied it to the asteroid- 
Jupiter system, again assuming zero eccentricity for 
Jupiter. Note that all these works used the secular ap- 
proximation for somewhat high a, with a\ — 3 AU; for 
this value the secular approximation breaks when the as- 
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teroid's apo-center crosse s Jupiter's orbit. In addition, 
IKinoshita fc Nakail (|2007l ) assumed a x = 1.814 AU, and 
as shown in the numerical integration in Figure [8] using 
the secular approximation for even this large value of a, 
may result in in-accurate evolutionary behavior. 

The Kozai-Lidov mechanism has been applied to the 
study of the outer solar system. Kinoshi ta fe Nakal 
(|2007|) also applied their analytical solution to Nep- 
tune's outer satellite Laomcdeia. This system has f^O 
and thus the TPQ li mit there is justified. In addition, 
iPerets fc Naoz! (|2009f ) have studied the evolution of bi- 
nary minor planets in the frame work of TPQ. In this 
problem e — > a nd thus the resul ts p resented there 
and la ter applied in lNaoz et al.l (|2010l) and lGrundv et al.l 
(120111) justify the applic ation of the TPQ limit. 

iLidov fc Ziglinl (J19761 sections 3-4) also solved ana- 
lytically the quadrapole-le vel approximation but, unlike 
IKinoshita fc Nakal (|2007| ). they did not restrict them- 
selves to the TPQ limit, and used the total angular mo- 
mentum conserva tion law in order t o calc ulate the in- 
clinations. Later, iMazeh fc Shahaml ([1979D also did not 
restrict their derivation to the TPQ limit (their eqs. Al- 
A8), and allowed for small eccentricities and inclinations 
o f the outer bod y. 

iKiseleva et al.l (119981) and 

lEggleton fc Kiseleva-Eggletonl (120011) st udied the 
Algol triple system ([Lestrade et al.l Il993l) using the 



TPQ equations. The TPQ equations were also used 
in the pape r that introduced the in f luential KCTF 
mech anism (jMazeh fc Shahaml 119791 : [Egglcto n et al.1 
Il998t ) . Figure Q2] compares the evolution computed in 
the (incorrect) "standard" quadrupole formalism, the 
correct quadrupole formalism, and the octupole- level 
forma lism applied to the Algol system of IKiseleva et al.l 
(|1998| ). The correct quadrupole formalism decreases 
the minimum value of 1 — ei by almost a factor of 2 
relative to the previous "standard" formalism. The 
reduced pericenter distance would strongly increase 
the effects of tidal friction (not included here), which 
may lead to rapid circularization of the inner orbit. 
The octupolc-lcvcl computation decreases the minimum 
p ericenter distance by a furth er 40%. 

iMikkola fc Tanikawal (|1998t ) analyzed the system CH 
Cygni, assuming that it was a triple svstenFH. They used 
the TPQ limit to model the system and derive its orbital 
parameters. They found that the best fitted model has 
in fact comparable masses for the inner and outer or- 
bits, and shorter periods for the outer and inner orbits 
then found in th e literature us i ng sp e ctroscopy and inter- 
ferom etry (e.g.. iHinkle et all 120 09; Mikolai cwska et all 
120101 ) . As we showed in Figure \W\ using the octapule- 
level equations we found a dramatically different evo- 
lutionary path for their best-fit parameters then the 
quadrapole-level, and also the TPQ. It may even be that 
th e triple model for this system cannot be excluded based 
on IMikkola fc Tanikawal (J19981 ) results, since the evolu- 
tion is so different under the corrected formalism. We 
s uggest that this wor k shou ld be r epeated. 

IWu fc Murravl (f2003h IWu et al.l (12007 



iFabrvckv fc Tremainel (|2007l ) and iCorreiaeTaTI ([20111 ) 
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11 It was re c ently claimed by Hinklc ct al] 1120091 ) ; 
IMikolai cwska ct al. (2010) that this system is in fact not a 
triple. 



Fig . 12. — The time evolution of the system Algol ( Egglcton ct al. 
1199811 . with (mi,m 2 ,m 3 ) = (2.5,2,1.7) M Q . The inner orbit has 
ai = 0.095 AU and the outer orbit has a 2 = 2.777 AU. The initial 
eccentricities are e± = 0.01 and e 2 = 0.23 and the initial relative 
inclination i = 100°. The z-components of the inner and outer 
orbital angular momentum, H± and ff 2 are normalized to the total 
angular momentum. The initial mutual inclination of 100° cor- 
responds to inner- and outer-orbit inclinations of 91.6° and 8.4°, 
respectively. We consider the (corrected) quadrupole-level evolu- 
tion (blue lines), octupolc-lcvcl evolution (dashed lines) and also 
the standard (incorrect) quadrupolc-lcvcl evolution. In the latter 
we have assumed, as in previous papers, that itot = ii, which 
results in the discrepancy between the inclination values. 



studied the evolution of a Jupiter-mass planet in stellar 
binaries in the framewo r k of KCTF. The case of HP 
80606 b (|Wu fc Murravl (I2003T); IFabrvckv fc Tremami 
(|20"07l their Fig. 1) and iCorreia et al.l (|201lL also Fig. 
1)) was considered with an outer stellar companion at 
1000 AU, and thus even if it that companion is assumed 
to be eccentric, e^/ is negligible, and the system is 
well described TPQ equations. Howev er, the sta t istica l 
distr ibution for closer stella r binar ies in IWu et"al] (|2007[ ) 
and IFabrvckv fc Tremainel (|2007l ) is only valid in the 
approximation where the outer orbit's eccentricity is 
zero. In fact for the eccentric and packed systems 
considered in those studies, cm is not negligible, and the 
the octupolc-lcvcl approximation results in dramatically 
different behavior (see §4.3). The same dramatic differ- 
ence i n behavior also exists in the a n alysis of triple stars 
(e.g., IFabrvckv fc Tremainel 120071 : IPerets fc Fabrvckvl 
120091 ) . and thus we suggest that these studies should be 
repeated 1 ^!. 

7. CONCLUSIONS 

We have shown that the "standard" Kozai formal- 
ism had an err or in the implementation of Hamilto- 
nian mechanics {Kozai 19621: ILidov! Il962h . Correcting 
the formalism we find that the ^-components of the both 
the inner and outer orbits' angular momenta in general 
change with time at both the quadrupole and octupole 
level. The conservation of the inner orbit's z-component 

12 We note that if em is too big the system becomes unstable, 
which may suggests that the phase space for which the eccentric 
Kozai is significant may be somewhat limited. 



Secular Dynamics in Three-Body Systems 



13 



of the angular momentum (the famous yl — e\ cos i = 
constant) only holds in the quadrupole-level test particle 
approximation. We have explained in details the source 
of the error in previous derivations f §6.1|) . 

We have re-derived the secular evolution equations for 
triple systems using Hamiltonian perturbation theory to 
the octupole-level of approximation (Section [5] and Ap- 
pendix El [4] and Appendix |C|) . We also showed that one 
can use t he simplified Ham iltonian found in the litera- 
ture (e.g.. iFord et al.ll2000bt ) as long as the equations of 
motion for the inclinations are calculated explicitly from 
the total angular momentum. 

The correction shown here has important implications 
to the evolution of triple systems. We discuses a few 
interesting implications in Section [5] We showed that 
already at the quadrupole-level approximation the ex- 
plicit assumption that the vertical angular momentum 
is constant can lead to erroneous results, see for exam- 
ple Figure 2] In this Figure we showed that far from the 
test particle limit in the quadrupole-level one can already 
find a significant difference in the evolutionary behavior. 
The corrected dynamics converges with the test parti- 
cle in the limit where GijG\ < 10 -4 , (see Figure [5]). 
During the evolution the inclination and eccentricity of 
both orbits oscillate. We show in Appendix IB1 that at the 
quadrupole level of approximation, the inner eccentricity 
and the mutual inclination have a well defined maximum 
and minimum. At the test particle limit these values 
converge to the critical inclination (39.2° > io < 140.8°) 
for large oscillatory amplitudes. 

We have derived the complete set of equations for the 
octupole-level evolution, including the explicit equations 
of motion for the evolution of the inclinations and the 
z-component of the angular momentum of the inner and 
outer orbits. 



The most notable outcome of the results presented 
here happens in the octupole-level of approximation, 
when the inner orbit flips from progradc to retrograde 
with respect to the total angular momentum (we call 
this flip th e "eccentric Kozai mech a nism" ) . We point 
out th at iKrvmolowski fc Mazehl (119991) ; IFord et al.1 
(I2000bfl: IBlaes et all (120021) : ILee fc PealeT (I2003bfl and 
lLaskar fc Bouel (|20100 had the correct equations of mo- 
tion, and could, in principle, have observed this phenom- 
ena. However, it seems that the assumption of a constant 
vertical angular momentum was built into the commu- 
nity understanding of Kozai mechanism that the eccen- 
tric Kozai-Lidov effect was overlooked. 

In lNaoz et al.l (|2011| ) we suggested that this effect may 
play an important role in the formation mechanism of 
retrograde Hot Jupiters. There we showed the impor- 
tance of this effect in a verity of planetary and stellar 
triple systems. Fo r some examples see Figures [5HS1 and 
iNaoz et all (|2011|) where we specifically discussed the 
evolution of two planet systems, triple stars and asteroids 
due to gravitational perturbations from Jupiter. We also 
compared our derivation with direct N-body integration 
and illustrated the same qualitative evolution. We also 
emphasized the importance of higher-orders approxima- 
tions, where cm is significant. 
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APPENDIX 
THE SECULAR AVERAGING AT THE QUADRUPOLE LEVEL 

Wc develop the complete quadrupole-level secular approximation in this section. As mentioned, the ma in difference 
between the derivation shown here and those of previous studies lies in the "eli mination of nodes" (e.g., Kozalll962t 
Ueffervs fc Moserlll966| ). which relates to the transition the invariable plane (e.g.. lMurrav fc Dermottll2000f ) coordinate 
system, where the total angular momentum lies along the z-axis. 

Transformation to the Invariable Plane 

We choose to work in a coordinate system where the total initial angular momentum of the sys tem lies along the z 
axis (see Figure [2]),; the x-y plane in this coordinate system is known as the invariable plane (e.g.. [Murray fc Dermott] 
l2000r i . and therefore we call this coordinate system the invariable coordinate system. We begin by expressing the 
vectors ri and r 2 each in a coordinate system where the pcriapse of the orbit is aligned with the x-axis and the orbit 
lies in the x-y plane, called the "orbital coordinate system," and then rotating each vector to the invariable coordinate 
system. The rotation that takes t he position vector in the o rbital coordinate system to the position in the invariable 
coordinate system is given by (see lMurrav fc Dermottll2000t chapter 2.8, and Figure 2.14 for more details) 



ri. 



R z {hi)R x (ii)R z {gi)ri, olb 



(Al) 



where the subscript "inv" and "orb" refer to the invariable and orbital coordinate systems, respectively. The rotation 
matrices R z and R x as a function of rotation angle, 9, are 



and 



(cos 9 — sin 9 
sin6l cos (9 
1 

(I ^ 

R x (9) = cosfl -sin6> 
V sin 9 cos 9 J 



Thus, the angle between ri and r 2 is given by: 

cos$ = rl mh R- 1 (g 2 )R~ 1 {i2)R z 1 (h2)R z {hi)R x {ii)R z {gi)r horh , 
where f 1,2.01b are unit vectors that point along ri^.orb- In the orbital coordinate system, we have 



I"l.2,orb 



'cos(Ji l2 )> 

sin(Zi j2 ) 

v , 



(A2) 



(A3) 



(A4) 



(A5) 
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Note that R~ x {h 2 )R z (hi) = R z {h\ — h-z) = R z (Ah), eso the Hamiltonian will depend on the difference in the longitudes 
of the ascending nodes; in a similar manner, the Hamiltonian depends on l\ and I2 only through expressions of the 
form li + gi and I2 +52- Replacing cos$ in the Hamiltonian, cq. (|T5j) . with equation (|A4[) we find a general form of 
the Hamiltonian in terms of the orbital elements. 

Transformation to Eliminate Mean Motions 

Because we are interested in the long-term dynamics of the triple system, we now describe the transformation that 
eliminates the short-period terms in the Hamiltonia n that depend of l\ and l 2 . The technique we will use is known as 
the Von Zeipel transformation (for more detalis, see lBrouwerlll959l ) . 

Write the triple-system Hamiltonian in eq. (fT5j) as 

H = H? + H 2 X +H 2 , (A6) 

where W'i and H^ are the Kepler Hamiltonians that describe the inner and outer elliptical orbits in the triple system 
and %2 describes the quadrupole interaction between the orbits. Note that H2 is O (a 2 ), and is the only term in H 
that depends on l\ or l 2 - Accordingly, we seek a canonical transformation that can eliminate the h and I2 terms from 
R.3. Such a transformation must be close to the identity, since H3 C?i; let the generating function be 

2 
S(L*,G*, Hj^g^hj) = ]T [Lft + Gtgj+Hjhi] +a 2 S 2 (L*,G*,H*,l j ,g j ,h j ), (A7) 

i=i 

where we indicate the new momenta with a superscript asterix, and S2 is the non-identity piece of the transformation 
that we will use to eliminate 'H.i . The relationship between the new and old canonical variables is 

OS „ 2 <9S2 ,. . 

Pi = «- = Pi + a -^— (A8) 

oqi aqi 

and 

* "& 9 0S2 , . „ % 

qt =wr qi+a w m 

where the momenta pi € {Li,Gi, Hi}, and the coordinates qi G {h,gi,hi}. Because our generating function is time- 
independent, the new and old Hamiltonians agree when evaluated at the corresponding points in phase space: 

nq i ,p i )=H*(q*,p*) (A10) 



when the phase space coordinates satisfy equations (|A8|) and (|A9|) . Inserting these relations into the un-transformed 
Hamiltonian, and expanding to lowest order in a , we have 



2 &HdS 2 2 &HdS 2 



"<*•*» + a dp^ -Qt - a d^ di = ^ M>M ■ (AU) 



and 



Equating terms order-by-order in a gives 

n?(q;, P ;) = Hi K (q;,p;), (A12) 

n^{qhp*)^n* 2 K {qhp*), (A13) 

i—l i—1 l 

Since the last two terms on the left-hand side of this latter equation are already O (a 2 ), only the "Hf and H^ parts 
of % contribute. These Kepler Hamiltonians only depend on h\ and L 2 , so there arc only two non-zero partials of % 
at order a 2 : 

nj 1 * *\ 1 2 9R-i dS 2 2 d% 2 9S 2 /Mt-N 

«»(*>*) + « ^^ +««£-«£- =**(«<.*)• (A15) 

We must use the terms that depend on S 2 to cancel any terms in H 2 that depend on I* and l 2 . Note that H2 is 
periodic in I* and l 2 with period 2ir (see equations (|A4[) and (|A5[l ). so we can write 

00 

n 2 (q*, P *)=a 2 h + a 2 J2 hk lk2 e- ikll *- ik2l K (A16) 

fel,fe2=l 



with 



ik lk2 = -^ f " dlldll n 2 (q*,p*) e *i'I+*fa«5. (A17) 
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Now let dH± jdL\ = wi(Li), and dH 2 /dL 2 = u^C^)- Suppose that S 2 is periodic in l\ and 1% (which arc equivalent, 
at lowest order, to l\ and 1%). Then 



a 2 h + a 2 ^2 h klk2 e- lkll i- lk2L i +a 2 to 1 ^ -ifcis fclfc2 e^ 



ik±li —ih<zli 



+ c?uj 2 Y, -ik2S klk2 e~ lk i h ~ lk * h = %l (q*,p*) , (A18) 

fel ,^'2 — 1 

where 

00 

S 2 = s + J2 s kl k 2 e- iklh - ik2h . (A19) 

fcl,fc2 = l 

The terms dependent on l\ will be eliminated from % 2 if 

s klk2 = -i , hk ]_ k2 , ■ (A20) 

Wlfci + w 2 fc 2 

Assuming than the system is far from resonance (that is, that U\k\ + u 2 k 2 7^ for all k\ and £2)1 this gives us the 
necessary S 2 to eliminate all terms in H 2 that depend on l\ or I2, leaving 

nl (q*,p*) = a 2 h = 4^2 / " ' dlldlin 2 (g*,pf)e iW *+*^. (A21) 

That is, our canonical transformation to eliminate the rapidly-oscillatin g p arts of T-L has left us with a Hamiltonian 
that is the average over the oscillation period of the original Hamiltoniarr^l. 
The value of the Hamiltonian in equation (|15|) averaged over the mean motions is 

^2 = ^rl [1 + 3 cos(2i 2 )] ([2 + 3e 2 ] [1 + 3 cos(2*i)] (A22) 

o 

+ 30e 2 cos(2 ffl )sin 2 (i 1 )) + 3cos(2A/i)[10e 2 cos(2pi) 
x (3 + cos(2ii)) + 4(2 + 3e 2 ) sin(ii) 2 ] sin 2 (i 2 ) 
+ 12(2 + 3e 2 - 5e 2 cos(2 3 i)) cos(A/i) sin(2zi) sin(2i 2 ) 
+ 120e 2 sin(ii) sin(2i2) sin(2<7i) sin(A/i) 
- 120e 2 cos(ii) sin 2 (i 2 ) sin(2gi) sin(2A/i)} , 
where C2 was defied in equation (|20|) . 

MAXIMUM ECCENTRICITY AND "KOZAI" ANGLES IN THE QUADRUPOLE APPROXIMATION 

First note that setting e'i = also means that G\ = 0. The values of the argument of periapsis that satisfy these 
relations are: g± = + wn/2, where n = 0, 1, 2... . Also, setting Gi(ei imaxm i n ) = means that i?i (ei ]maXjln i n ) = and 
i\ = 0, i.e., an extremum of the eccentricity is also an extremum of both the inner and outer inclinations. 

The conservation of the total angular momentum, i.e., Gi + G2 = G to t sets the relation between the total inclination 
and inner orbit eccentricity. We re-write equation ([6]) as 



L{{1 - e{) + 2L X L 2 yjl -e 2 ^/l-e 2 cosi tot = G 2 tot - G 2 2 , (Bl) 

where in the quadrupolc- level approximation e 2 and G 2 are constant. The right hand side of the above equation is set 
by the initial conditions. In addition, L\, and L 2 [see eqs. ((3]) and (j4])] are also set by the initial conditions. Using the 
conservation of energy we can write, for the minimum eccentricity case (i.e., setting g\ = 0) 

-^ = 3cos 2 i tot (l-e 2 )-l + 6e 2 , (B2) 

where we also used the relation A/i = w. We find a similar equation if we set g\ = tt/2: 

-^ = 3cos 2 i tot (l + 4e 2 )-l-9e 2 . (B3) 

ZC 2 

13 Note that the canonical variables are also transformed. They , .. , ., , , n , ., . . . . ■ , , m i 2\ j 

,.„ . ,, ■ ■ 1 ■ , 1 . /n r i\ tt , , ■ j-c orbits described by fii, as this interaction is already O \a , and 

diner from the original variables at Ota ). However, this dil- ,, jrr v *. 4 u ■ ■ 1 j 4. c j v '• ui 

e ... f , , . . . . v . / , ' . , so the differences between the original and transformed variables 

tcrcncc is irrelevant when evaluating the interaction between the . ., , . , , ,. , 

& contribute at sub-leading order. 
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Fig. 13. — The total inclination and eccentricity rela tion . We show constant energy curves (solid curves, the "ha lf" e llipses are for eq. | |B5| | 
and two examples for io = 80° and io = 100° of eq. I IB6I 0. and constant total angular momentum curves (eq. JB4t dashed curves). The 
initial conditions considered here are ej = 0, gj an d e 2 = " arlc ' L1/L2 = 0.07, appropriate for the Algol system (see £|5.4| I, We consider 
four different initial inclinations and their symmetric 90° counterparts, from bottom to top 10, 30, 60 and 80 degrees. We also show an 
example (highlighted curve) for the Algol system which is a result of integration of the quadrupole-levcl approximation equations. 



Equations (|B1|) . (|B2[) and (|B3[) give a simple relation between the total inclination and the inner eccentric ity. The 
remainder of the parameters in the equations are defined by the initial conditions. Thus, using equatio ns ()B2|) and 
(|B1[) we can find the minimum eccentricity reached during the oscillation and using equations (|B3[) and (|B1|) we can 
find also the maximum and the minimum inclinations. The following example illustrates the relation defined by these 
equations between the inclination and the eccentricity. 

For simplicity we set initially e\ = 0, g® and e® = (the superscript stand for initial values). In this appendix we 
consider only the quadrupole-level approximation, and thus e 2 doesn't change. Using these initial conditions (and for 
some initial mutual inclination ?o) we can write equation (|B1|) as 



1 — e 2 cos it t = cos io 



_£i_ 

2L 2 



(B4) 



We show these curves for different i$ in Figure [TBI ( short dashed curves) for a hypothetical system with the parameters 
of Algol (but with e 2 = 0, see £|5.4[) . Note that there is a slight asymmetry b etween the prograde and retrograde orbit s 
due to the L\/L 2 factor (whi ch is not t he case for the test particle case, see[Li thw ick fc Naodl201lt iKatz et al.ll2011l ). 
We also write equations (|B2|) and (|B3| using the initial conditions. Equation (|B2j) can be simplified to 



(1 



1 cos itot = cos i + 2e\ 



(B5) 



depicted in Figure [T3] (solid curves, for different io). As ca n be seen from the Figure, this equation gives the minimum 
eccentricity, which is the crossing point with equation (|F34|) . For these choice of initial conditions the minimum 



eccentricity is e" = 0. Equation (|B3j) becomes 



(4 



-t-ejeos i to t 



cos 2 io 



3e? 



(B6) 



whic h is depicted in Figure [13] (long dashed curves, for io = 80° and 100°). We now use this equation and equation 
(JB4I) to find the maximum eccentricity. After some algebra we find: 



3 + 4-— cos i 
L 2 




2L 2 



-3 + 5 cos 2 i = 



(B7) 
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As we approach the TPQ limit, L2 ^> L±, and this equation becomes 

e\ = 1 - - cos 2 i , (B8) 

which gives the maximum eccentricity as a function of mutual initial inclination with zero initial inner eccentricity. 
In Figure [131 w e show that this approximation still hold fairly well even for the Algol system, where L\jL 2 ~ 0.07. 
Equa tion (|B8|) has been found previously (e.g. Ilnnanen et al. I [19971 iKinoshita fc Nakal 119991 : iValtonen fc Karttunenl 
12006( 1 in the TPQ approximation, but in these works it is assumed to valid outside that limit. A solution exists only if 
the right hand side of this equations is positive, thus we find the critical angles for large Kozai oscillation in the TPQ 
limit: 

39.2° > i < 140.8° . (B9) 

For larg er L\jL 2 and/or for initial e\ > this limit and e max arc different and the full solution of equations (|B1[) JF32|) 
and lB3p is required. In fact for each initial set of e\ > and i toi , there is a specific Li/L 2 that will produce an angular 
momentum curve that crosses 90°. Thus, for initial g\ > 90° the mutual inclination can oscillate from value below 90° 
to above. This happens because the inclination of the outer orbit i 2 changes considerably, while the inner orbit remain 
prograde (if started prograde). We emphasize that the "ocatpolc-Kozai" behavior (2J) is of course present; and can 
only be neglected when cm <C 1. 

THE FULL OCTUPOLE-ORDER EQUATIONS OF MOTION 
We define: 

15 fc 4 (mi + m 2 ) 9 m\{mi-m 2 ) L\ 
3 ~~16T(mi+TO2+TO 3 ) 4 (mim 2 ) 5 L\G\ ' [ ' 

Note that t his d efinition is with a different sign from iFord et al.l (|2000bh . and consistent with iBlaes et al.l ((2002); 
iFord et aLl (j2004f ). For equal mass m\ and 777.2 this factor is zero. We also define: 

^ = 4 + 3e 2 -^sinz t 2 ot , (C2) 

where 

B = 2 + 5e 2 -7e 2 cos(2 3 i) , (C3) 

and 

cos <j) = — cos (71 cos g-2. — cos z to t sin g\ sin g 2 ■ (C4) 

The time evolution of the argument of periapse for the inner and outer orbits are given by: 

.91 =6C 2 |^-[4cos 2 i tot + (5cos(2 gi ) - 1) (C5) 

x (1 - e 2 - cos 2 H ot )] + £2^2l[ 2 + e \{Z - 5 cos(2 9l ))] 

O2 

C 3 e 2 { ei ' 



G 2 G\ 

x [sin 5l sin 5 2 (10(3 cos i 2 ot - 1)(1 - e 2 ) + A) 

1 - e 2 

- 5 cos itot cos </>] -p- x [sin gi sin g 2 

eiGi 

x 10 cos itot sin i 2 ot (1 - 3e 2 ) 
+ cos (j)(3A - 10 cos i 2 ot + 2)] 



and 



g 2 = 3C 2 



1 
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(C6) 



■£- [4 + 6e 2 + (5 cos 2 i tot - 3) (2 + e 2 [3 - 5 cos(2 5l )] ) 
f4e 2 _+l. 

G 2 



-C 3 ei<j singisin^l " 2 „ 1 10cosi t ot sin 2 H ot (l - ef) 



e 2 



-cost 



Gi 



cos l tot 



[A+ 10(3 cos 2 itot-l)(l-e 2 )] 



5i?cosi to te2 ( -^~ 



COSltot 

G 2 



4e|±i A 

e 2 G 2 



The time evolution of the longitude of ascending nodes is given by: 

3G 2 



hi=- 



(2 + 3e 2 - 5e 2 cos (2gi)) sin (2z to t) 



G\ sinii 
— G 3 eie 2 [5-Bcosi to tcos0 — A sin 51 sin 32 + 10(1 + 3cos 2 i to t) 

X (1 — ei) singi sing 2 ] 



(C7) 



Gisinii 



where in the last part we have used again the law of sines for which sinii = G 2 sini tot /G t ot- The evolution 
longitude of ascending nodes for the outer orbit can be easily obtained using: 



The evolution of the eccentricities is: 



ei = G 2 



Gi 



hi 



i[30ei sin 2 itot sin(2gi)] 



of the 

(C8) 

(C9) 



+ G 3 e 2 1 [35 cos </> sin 2 i to te 2 sin(2g-i) 

Gi 

-10 cos itot sin 2 i to t cos g\ sin g 2 (l - e 2 ) 

- A(sin 5! cos g 2 - cos i tot cos g x sin g> 2 )] , 



and 



<?2 : 



-G 3ei - 



G 2 



■[10 cos itot sin 2 i tot (1 - e 2 ) singi cosji 2 



- A(sin gi cos g 2 - cos i to t cos gi sin g 2 )] ■ 
We also write the angular momenta derivatives as a function of time; for the inner orbit 
Gi = -G 2 30ejsin(2.g 1 )sin 2 (i t ot) + G 3 eie 2 ( 

— 35e 2 sin 2 (z to t) sin(2<7i) coscj) + A[sin(gi) cos(g 2 ) 

-cos(i to t)cos(gi)sin(g 2 )] + 10cos(i to t) sin(i to t)[l - e 2 ] cos(fl-i) sin(g 2 )) , 
and for the outer orbit (where the quadrupole term is zero) 

G 2 = C 3 eie 2 [A{cos(gi) sin(g 2 ) - cos(i tot ) sin(gi) cos(g 2 )} 
+ lOcos(itot) sin 2 (i to t)[l - e 2 ] sin(g>i) cos(g 2 )] . 
Also, 



(CIO) 



(Cll) 



tt - Gl r G2 r 

tl\ — — (-T1 — — (_T 2 



where using the law of sines we write: 
The inclinations evolve according to 



Hi = - 



Gtot 

smi2 



Gi 



• srrui • 

l^l : — : — ( - r 2 

sm itot 



sin itot 



r ' • \ Hi Gi 

(cosii) = — — - COSli , 



(C12) 

(C13) 
(C14) 

(C15) 
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and 

(cosi 2 ) = -^ - -^- cosi 2 . (C16) 

(_r 2 2 

Our equations are equivalent to those of iFord et al.l (|2000bD . but we give the evolution equations for Hi and il 2 (and 
i\ and i 2 ). 

ELIMINATION OF THE NODES EQUATIONS DERIVATION 

Here we show a derivation of equations (|D7[) which lead to the error in previous treatments. The angular momentum 
vectors of the two orbits are given by 

Gi j2 = Gi,2 (sin ii ]2 sin fti^, — sin 11,2 cos /11,2,00s ii )2 ) ■ (Dl) 

Thus the total angular momentum vector is then: 

Gtot = [G\ sinii sin/ii + G2 sini 2 sinft, 2 , (D2) 

— G\ sinii cos/ii — G 2 sinz 2 cos/i 2 ,Gi cosii + G 2 cosi 2 ) . 

Recall that the z-component of the angular momentum is Hj = Gj cos ij . 
In the elimination of the nodes we set Gtot||-z thus, 

Gi sinii sin/ii = — G 2 sini 2 sin/i 2 , (D3) 

Gi sinii cos/ii = — G 2 sini 2 cos/i 2 , (D4) 

H 1 +H 2 = G tot . (D5) 



Dividing Eqs (|D3|) and (|D4|) . we obtain 
implying that 



tan/ii = tanft. 2 , (D6) 

hi - h 2 = 7r . (D7) 

Because the dynamics of the system conserves total angular momentum, this result will always hold. This is, however, 
a dynamical restriction, and does not imply any restriction on the partial derivatives that produce the equations of 
motion from the Hamiltonian. In other words, we cannot substitute 

hi = h 2 - 7T . (D8) 

into the Hamiltonian before computing equations of motion that may involve terms with 

£■ 

EXAMPLE 

As a simple example to illustrate the source of the mistake, let us consider a 1-D system with two equal masses 
connected by a spring. The Hamiltonian for such a problem is 

H = &&l k -^-^- < E1 » 

where k s is the spring constant. Qualitatively equivalent to the elimination of the nodes would be here to transform 
to the center of mass of the system, so that xi = — x 2 . If we now substitute this relationship between the coordinates 
into the Hamiltonian, we get 

But this is incorrect! This Hamiltonian implies, for example, that P 2 = const. 

Note that the error that leads to the incorrect secular three-body Hamiltonian is analagous: conservation of mo- 
mentum gives a relation between two coordinates (hi = /i 2 — n), and substitution of this relation into the Hamiltonian 
gives the incorrect relation ^J\ — ei cosii = const. 



